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ABSTRACT 

A  three  dimensional  primitive  equation  model  based  on 
the  Boussinesq  equations  is  developed  and  applied  to  the 
mesoscale.   The  model  is  tested  with  one  dimensional  flow  in 
a  balanced  state  and  with  two  cases  of  gravity  waves,  which 
are  forecast  for  30  minutes  and  compared  with  the  analytic 
solutions. 
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I.   INTRODUCTION 

Several  recent  attempts  to  understand  the  complex 
dynamic  interaction  of  the  thunderstorm  and  its  environment 
have  been  made  using  the  dynamical  thunderstorm  model  pre- 
sented by  Newton  and  Newton  (1959) .   The  model  proposes  that 
strong  vertical  drafts  within  the  thunderstorm  create  an 
effective  barrier  to  the  environmental  flow.   The  analysis 
by  Browning  (1964)  concerning  airflow  near  and  within  severe 
thunderstorms  provides  valuable  additional  knowledge  about 
the  character  of  the  vertical  motion.   Convective  storm  mass 
and  water  budget  assessments  such  as  the  one  conducted  by 
Newton  (1966)  substantiate  the  concept  of  upper  level 
entrainment  of  cool,  dry  air  first  recognized  by  Normand 
(1946)  . 

Radar  determined  chaff  trajectories  were  used  success- 
fully by  Fankhauser  (1968)  to  follow  the  motion  of  air 
parcels  in  and  around  a  thunderstorm.   Further  use  of  chaff 
will  undoubtedly  make  possible  a  good  understanding  of  the 
paths  taken  by  environmental  parcels  from  all  positions 
around  a  local  storm. 

As  noted  by  Lilly  (1962) ,  the  evolution  of  convective 
motions  with  relatively  large  amplitudes  can  best  be 
described  through  the  combined  application  of  conventional 
analytic  methods  and  numerical  experimentation.   Ogura  (1963) 
discussed  the  possible  use  of  the  primitive  equations  as  a 
basis  for  forecasting  small  scale  dynamical  features.   He 


11 


suggested  that  interaction  between  convective  elements  and 
the  prevailing  wind  field  could  thus  be  investigated  numeri- 
cally.  The  results  of  such  efforts,  correlated  with  the 
subjective  modeling,  should  produce  verifying  as  well  as 
augmenting  information  for  the  subjective  products. 

A  numerical  study  of  turbulent  flow  conducted  by 
Deardorff  (1970)  makes  use  of  a  three  dimensional  primitive 
equation  model.   The  forecast  procedure  he  uses  is  similar 
to  the  procedure  used  in  this  study. 

The  numerical  model  presented  here  is  a  three  dimen- 
sional primitive  equation  model  based  on  the  Boussinesq 
equations.   Forecasts  are  made  for  the  wind  components  and 
potential  temperature  after  imposing  initial  conditions  on 
these  same  parameters. 

The  grid  used  for  computations  is  25  by  25  in  the 
horizontal,  with  15  levels  in  the  vertical.   Grid  interval 
in  all  directions  is  1  kilometer.   A  mean  latitude  of  35° 
is  used  to  determine  the  Coriolis  parameter,  which  is 
treated  as  a  constant. 

Predictions  are  made  using  the  IBM  OS/360  computer  of 
the  W.  R.  Church  computer  center.   A  one  minute  time  step 
was  chosen  to  conform  to  computational  stability  require- 
ments. 
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II.   DEVELOPMENT  OF  THE  MODEL 

A.   THE  FORECAST  EQUATIONS 

The  pressure  term  in  the  equation  of  motion  may  be 
expressed  as 


1   «      RT  „ 
P       V3  P  =  T    V3 


5        1/V 


i/fc 


*>V3j3  , 


where 


and 


&    =  ct 


1  k 


o 


0  =   T 


o 
P 


Using  this  result  in  the  component  equations  of  motion 


yields 


du 
dt 


=  -  © 


a/3 


X 


+  fv  +   ]/V3   u, 


(1) 


dv 
dt 


=  -  © 


a/3 


-  fu  +   VV3   v, 


(2) 


dw 
dt 


=  _  0  AfL    _  g  +     yv2    W. 


(3) 


Forming  the  perturbation  expressions 


/S  =  /50(z)  +  /3'  (x,y,z,t) 


and 
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e    =       $0  +     0  (x,y,z,t) 
where 

[ £'|    «   /5o    and  |  o  \«e0  ' 

and  substituting  into  (1)  yields 

#  =  -  0O  -M*  +  fv  +  V  V3   u  ,  (4) 

ox 

ofi' 

where  the  relatively  small  product  0  rr  —  ■   is  neglected. 

o  x 


Similarly, 


dv    .  *   _§JL  _  f u  +  v  Vo2  v.  (5) 


—  —  ft  _  -r       —  x  u  -r   f  v  o 

dt     "°  dv  3 


For  the  vertical  component  equation  the  substitution  yields 


&  =-*o!f0  -  •o-H"  -•■£  -9+VV32v, 


Since    /S         and       0  are    related  by    the  hydrostatic    equation 

oflo     =   _  JL 
d  z  #o 

the  vertical  component  equation  may  be  written 

*Z=   -  0  -Ml  +  i?  +   W  *  w.  (G) 

dt         *o  az        0O  3 

Defining    <t>    =    $     jS '   an<^  substituting  into  (4),  (5)  and 
(6)  results  in  the  system  of  equations 


14 


—     =      -£(u)     -  — —     +    fv  +      W,      u  (7) 

dt  ox  3  x 


Al     =      -£(v)     -  -|^      -    fu  +      W,      v  (8) 

at  ^y  3 


^    .    _£(w)  _A£   +  2i_   +    vv  2  w  (9) 


ot  dz  fl  3 


-f£     =      -£    (0)  (10) 

dt 


v3   •     v3     =0  (11) 


where   £ ( S)  =  -v3 ■  V3  5,8=  u,v,w,  0.   Equation  (11)  is  the 
continuity  equation.   In  order  to  eliminate  sound  waves, 
incompressible  flow  is  assumed.   Ogura  and  Phillips  (1962) 
have  shown  that  this  system  of  equations,  which  they  call 
the  anelastic  equations,  is  quite  accurate  for  shallow  con- 
vection. 

B.   DETERMINATION  OF  THE  PRESSURE  TERMS 

Equations  (7) ,  (8) ,  and  (9)  may  be  combined  and  written 
in  the  vector  form 

dVo  _,  _>  ^^90^  2  _ 

Y£     =      -(V3-  V3)V3    -  V3<t>  -f  (k  x  v3)    +  —  k  +  yv3     V3    . 

(12) 

Taking  the  three  dimensional  divergence  of  (12)  yields 

15 


(V3    •    V3) 

*  =   0 


at 

—      —    T1 


[(v3  •  v3)  v3J  _    v3  <p 


-V3    •       [f(Xx   V3)]+^- 


(13) 

be 


For   a    limited   horizontal    scale    the   Coriolis   parameter   may  be 
treated   as   constant.       Thus,     (13)    becomes 


V32<^     =      V3    '       [(^3    *    V3^3]~    fo   V3    '     (*  X  ^3)    + 


_9_  _ze 
6Q    az 


=  -  v. 


[<V3  •   V3)  V3]  +  fQ  C   +  £-§£   .   (14) 


which  can  be  solved  using  relaxation  techniques 
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III.   FINITE  DIFFERENCING 

Equations  (7)  through  (10)  are  solved  numerically  by 
introducing  finite  differences  in  x,  y,  z  and  t.   Solution 
for  the  pressure  term  can  be  carried  out  by  applying  the 
extrapolated  Liebmann  relaxation  technique  to  evaluate  (14) . 
Appendix  A  shows  the  development  of  the  relaxation  equation. 
The  procedure  used  to  compute  the  forcing  function  is 
included  in  Appendix  B. 

In  order  to  maintain  finite  differencing  consistency 
throughout  the  entire  forecast  process,  it  is  necessary  to 
use  the  finite  difference  forms  of  equations  (7),  (8)  and 
(9)  to  form  the  finite  difference  equivalent  of  equation 
(14). 

The  vertical  domain  consists  of  fifteen  levels  with 
values  of  u,  v,  w,  $,    and  <p    assigned  to  each  grid  point. 
The  space  differencing  method  used  is  based  on  a  scheme 
devised  by  Arakawa  (1966)  and  is  designed  to  conserve  total 
energy.   The  non-linear  operator,  £(s) ,  takes  the  finite 
difference  form 

£(s)    ■    [(ui+i,j,k+  ui.j,x)  (8i+i.j.*+  sij.k' 

"  {ui-i,j,*+  ui,j.k'    (si-i,j,x+si,j,k>]/4Ax 

+    [(Vi,J  +  l,*+    ^.j.k'     (Si.j  +  l,k+    Si.j,k» 

-    (Vi,j-l,k+    vijk)     (Sijj.lik+   Sijk)]     /4Ay 
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+   [<wi,j,k+l +  wi,j.k>  (si.j.k+l +  si,j,k) 

-  <wi,j.k-l  +  wi.j,k>  <si.j,k-l  +  si.j,k>]   /4A: 


The  linear  space  derivatives  are  of  the  centered  form 
£3      (Si+l,j,k  "  si-l,J,k> 


dx  2  Ax 

Centered    time   differencing    is    used    for   all    forecasts    except 
the    first  which   employs    a    forward    time    step. 
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IV.   BOUNDARY  CONDITIONS 

The  boundary  conditions  are  selected  to  enforce  mass  and 
energy  conservation  within  the  domain  of  the  grid.   Figure  1 
illustrates  the  grid  domain  and  boundary  placement  in  the 
horizontal.   Boundaries  in  the  vertical  are  similarly 
placed,  with  bottom  and  top  boundaries  at 

2+3 
fc  = =  2.5 


and 


(KM-1)  +  (KM-2) 


respectively,  where   KM   is  the  number  of  grid  points  in  the 
vertical. 


i=l 
j=JM 


*    « 

r— 


i=IM 
j=JM 


North  boundary 


West i boundary 


Eas  t  J  boundary 


4 1 


South  boundary 


i=l 
j=l 


Figure  1:   Horizontal  grid 


i=IM 


Periodicity  in  the  east-west  direction  ensures  cyclic 
continuity  of  all  parameters.   This  condition  is  imposed  by 
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applying  the  relationships 


Sl,j,k    ~  SIM-3/j/k  ' 


S2,j,k      SIM-2,j,k  ' 


IM-l,j,k  '   3,j,k 


IM,  j  ,k       4,  j  ,  k 


where   S  =  u,  v,  w,  0  and  <£.   Since  values  of   S   at  points 
in  the  outer  columns  are  required  in  solving  for  the  pressure 
term  and  in  calculating  divergence,  a  four  point  periodicity 
constraint  is  necessary. 

In  order  to  prevent  the  flow  of  mass  and  energy  through 
the  northern  and  southern  boundaries,  a  condition  of  no  flux 
is  imposed  on  the   v  component  of  the  wind.   This  is  accom- 
plished by 


v.  o  i   -   -v.  o  n 
i, 2,k       l, 3,k 


and 


Vi,JM-l,k      vi,JM-2,k   * 


The  values  of   v  on  the  outer  row  are  obtained  by 

V-i  i  v  =   v.  .  n 
1/  J-/K     i,4,k 


and  v .      =   v . 

i,JM,K     i,JM-3,k 
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The   u  and  w  wind  components  and  0      are  projected 
outward  by 


i, l,k  "   i,2,k  '   i, 3,k 


and  si,JM,k  "  Si/JM-l,k  "  Si,JM-2,k   ' 


where  s  =  u,w,  0  . 

Values  of  pressure  on  the  outer  rows  are  obtained  by  using 
the  geostrophic  relationship 

30  =      -fu 

Sy 

The  upper  and  lower  boundaries  are  treated  in  the  same 
manner  as  the  north  and  south  boundaries,  with  the  no-flux 
condition  being  imposed  on  w  by 

w.  .  0     =  -w 


±,3,2  -      wi,J,3   * 


Wi,  j,KM-l  "  ~Wi,j,KM-2  , 


w.  .  ,     =  w.  .  .   , 


and  w.        —  Wj  j  VM  o 

i,j,KM      i,j,KM-i 

Values    of      u      and      v     on   the    top   and  bottom   levels    are 
obtained  by 

Si,j/1   =    Si,j,2    =   Si,j,3 
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and 


Si,j,KM  "  Si,j/KM-1  :  Si,j,KM-2   ' 


where   s  =  u,  v. 

The  hydrostatic  relationship, 

d  z     0o 

is  used  to  determine  values  of  0    and  <j>     above  the  top  and 
below  the  bottom  boundaries.   The  application  of  proper 
boundary  conditions  was  found  to  be  critical.   Appendices  C 
and  D  show  the  development  of  these  conditions  in  detail. 
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V.   INITIAL  CONDITIONS 

Initialization  of  the   u,v,w   and  Q      fields  was 
accomplished  according  to  the  experiment  being  prepared. 
During  the  check-out  phase  of  programming  the  model,  a 
steady-state,  geostrophic,  zonal  flow  was  used  to  provide  a 
check  on  computational  stability.   Subsequent  tests  were 
initialized  to  include  gravity  waves,  both  stable  and 
unstable. 

A.   ZONAL  FLOW  CASE 

Stability  of  the  balanced  state  was  tested  by  inserting 
a  10  m  sec"   wind  throughout  the   u   field.   Initial   v 
and  w   fields  were  set  to  zero.   Standard  atmospheric 
values  of  pressure  were  introduced  and,  by  using  an  adia- 
batic  lapse  rate  of  temperature,  the  equation 


$     =      T 


po 


provided  the  initial  0     field.   Since  the  relaxation 
technique  used  to  obtain  the  initial  <f>     field  proved  to  be 
extremely  efficient,  choosing  the  initial  guess  <t>     field  as 
an  arbitrary  constant  was  convenient. 

B.   GRAVITY  WAVE  CASES 

Linearized  gravity  wave  solutions,  accurate  for  small 

w 

scale  perturbations  only,  were  used  to  test  the  model.   For 
the  stable  gravity  wave  the  following  relationships  describe 
the  initial  conditions. 
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u     =  — 


A  m  W  \ 

-o r     sin  Ax  cos    n  y      cos   mz 

A2  +  n2  " 


v       =  — 


fj  m  W 

-5 —     cos     Ax   sin    IJy      cos   itiz 

a2  +  \x2 


w 


=     W  cos  Ax  cos    [X  y      sin 


mz 


*  VI         b  9  •       \ 

Q     =  .—-     -^-^     sin  Ax  cos    i/y      sm  mz 

Ac   d  z 


Phase  speed  may  be  calculated  from  the  relationship 


N 
c   =  — 


A2  +  /i2 

A  *  A2  +  /i2  +  m2 


(15) 


For  the  gravity  wave  associated  with  unstable  stratifi- 
cation, the  initializing  relationships  are 


u     = 


v     = 


m 


m 


1  + 


1  + 


N 


n' 


N' 


n' 


W   sin  Ax  cos  fA  Y      cos   mz 


W  cos  Ax   sin  JJy    cos   mz 


w 


=      W  cos  Ax  cos   Uy     sin 


mz 


0  =~¥L     ..?" .     cos  Ax  cos  fjiy     sin  mz    , 


n       d 


where 


N2  =|-    -M 

^O  ^Z 
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and 

„2   _    f*\+i*    ,  as) 

where   n   is  the  growth  rate  of  the  disturbance. 

The   </>  field  may  be  initialized  for  the  gravity  waves 
exactly  as  in  the  zonal  flow  case. 
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VI.   RESULTS 

The  balanced  state  introduced  for  the  one  dimensional 
flow  case  as  a  check  on  the  computational  stability  and 
forecast  capability  of  the  model  should  remain  indefinitely, 
subject  to  the  validity  of  the  approximations  in  the  basic 
equations,  round-off  error  in  the  computer,  accuracy  of  the 
assumptions  at  the  boundaries,  and  stability  requirements. 
The  use  of  proper  boundary  conditions  was  the  most  important 
factor  in  determining  the  length  of  the  forecast  produced. 
Until  the  combination  of  boundary  conditions  shown  in 
section  IV  and  Appendix  C  was  utilized,  the  forecasts  were 
valid  for  no  longer  than  5  minutes.   Longest  forecasts  of 
the  balanced  flow  were  obtained  by  applying  the  boundary 
conditions  that  are  periodic  in   x,  geostrophic  in   y,  and 
hydrostatic  in   z. 

Accuracy  of  forecasts  was  monitored  through  the  calcula- 
tion of  the  three  dimensional  divergence,  which  should  have 
remained  very  small  throughout  the  prognosis.   The  smallest 
possible  relaxation  tolerance  consistent  with  single  preci- 
sion mode  capability  of  the  computer  was  found  to  allow  the 
least  amount  of  divergence.   When  larger  tolerances  were 
allowed,  a  substantial  amount  of  divergence  developed  on  the 
north  and  south  boundaries  at  mid-levels  and  increased  with 
height. 

A  test  was  conducted  to  determine  the  optimum  relaxation 
coefficient  for  the  grid  interval  and  grid  domain  which  the 
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model  employs.   The  value  selected  did  not  appear  too 
critical,  as  long  as  it  was  between  1.5  and  1.7. 

For  the  first  experiment  that  produced  a  relatively  long 
prognosis,  an  adiabatic  lapse  rate  of  temperature  and 
standard  atmospheric  values  of  pressure  were  used  to  ini- 
tialize the  potential  temperature  field.   This  produced  a 
slightly  unstable  average  lapse  rate  of  potential  tempera- 
ture, with  two  thin  stable  layers  included  near  the  bottom. 
With  these  conditions  the  model  forecast  for  47  minutes. 

The  time  step  was  then  decreased  from  1  minute  to  40 
seconds,  with  the  result  that  the  forecast  proceeded  for  72 
time  steps,  which  was  essentially  the  same  length  as  the 
previous  forecast.   In  an  effort  to  improve  the  initial 
balance,  pressure  values  which  would  lead  to  a  nearly  neutral 
lapse  rate  of  potential  temperature  were  inserted  in  the 
first  experiment.   The  forecast  continued  for  56  minutes. 
Returning  again  to  the  conditions  used  in  experiment  one, 
the  lapse  rate  of  potential  temperature  was  modified  to 
include  stable  layers  across  the  top  and  bottom  boundaries. 
With  these  conditions  the  model  produced  a  forecast  for  69 
minutes. 

Introduction  of  a  small  perturbation  in  the  form  of  a 
gravity  wave  provided  another  method  of  observing  the  model 
behavior.   The  gravity  wave  in  stable  stratification  should 
propagate  in  the   x  direction  only,  without  growing.   With 
unstable  stratification  the  wave  should  not  propagate,  and 
should  grow  at  a  rate  given  by  equation  (16) . 
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The  wave  inserted  for  the  stable  case  was  forecast  for 
30  minutes,  after  which  a  comparison  was  made  between  the 
initial  and  forecast  values  of  u,  v,  w  and  0   .   With  this 
comparison  it  was  possible  to  observe  the  amount  of  spurious 
growth  and  the  speed  of  propagation.   The  relationship  of 
the  initial  and  forecast  wave  forms  in   x   for  the  line  j=8, 
k=6  is  shown  in  Figures  2a,  2b,  2c  and  2d.   The  exact  solu- 
tion of  equation  (15)  indicated  that  the  phase  speed  for  the 
case  shown  was  10.2  m  sec   .   Phase  speeds  calculated  from 
the  plotted  values  were  10.5  m  sec   ,   10.4  m  sec-1  , 
9.5  m  sec    and  10.3  m  sec"    for  u,  vr  w  and  $     respec- 
tively.  The  difference  in  phase  speeds  between  w   and  0 
was  in  no  doubt  related  to  the  fact  that  the  wave  grew.   The 
wave  forms  in   y   illustrated  in  Figures  3a,  3b,  3c  and  3d 
show  that  the  wave  grew  in  the   y   direction  but  did  not 
propagate.   Figures  4a,  4b,  4c  and  4d  show  that  no  propaga- 
tion of  the  wave  occurred  in  the  z      direction,  but  that  the 
growth  was  not  smooth. 

For  each  combination  of  initial  conditions  applied  to 
the  wave,  growth  was  apparent.   By  modifying  the  lapse  rate 
of  potential  temperature  toward  adiabatic  conditions,  the 
growth  decreased. 

The  gravity  wave  with  unstable  stratification  was 
compared  in  the  same  manner  as  was  the  stable  case.   Figures 
5a,  5b,  5c  and  5d  show  that  the  wave  did  not  propagate. 
Furthermore,  the  growth  that  was  experienced  by  the  wave  was 
very  close  to  that  which  was  expected.   The  amplification 
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factor  based  on  the  computed  growth  rate  of  .1376  x  10~^ 
sec   was  determined  to  be  11.9  at  30  minutes,  and  amplifi- 
cation factors  of  11.8,  12o0,  12.4  and  11.3  were  noted  for 
u,  v,  w  and  $    respectively. 
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VII.   SUMMARY  AND  CONCLUSIONS 

A  multi-layer  primitive  equation  model  was  applied  to 
a  mesoscale  grid  for  two  test  situations.   First,  a  one 
dimensional,  balanced,  steady  state  current  was  initialized 
to  check  the  computational  stability  and  finite  difference 
form  of  the  equations  of  motion,,   Then  an  initial  gravity 
wave  perturbation  was  introduced  for  both  stable  and 
unstable  stratification. 

The  stable  stratification  case  yielded  a  reasonably 
accurate  phase  speed  of  propagation  but  some  undesirable 
growth  of  the  disturbance  was  observed,  as  was  discussed  in 
the  previous  section. 

The  unstable  stratification  case  produced  nearly  correct 
growth  rates  but  some  difference  in  the   w   and  0     growth 
rates  may  be  an  indication  that  a  problem  exists  similar  to 
that  in  the  stable  case,  where  w  and  0    waves  were  propa- 
gated at  slightly  different  speeds 0   Apparently,  some  small 
error  still  exists  in  the  program. 

The  fact  that  the  growth  of  the  propagating  gravity 
wave  can  be  decreased  by  inserting  a  more  nearly  adiabatic 
lapse  rate  of  potential  temperature  can  be  explained  by 
observing  that  the  buoyancy  term  in  the  Boussinesq  equations 
becomes  increasingly  inaccurate  with  increased  height.   By 
assigning  to  $       a  value  of  potential  temperature  from  a 
middle  level,  the  error  introduced  in  the  buoyancy  term  is 
minimized. 
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Since  the  Boussinesq  equations  are  accurate  only  for  a 
shallow  atmosphere ,  the  model  could  be  expected  to  forecast 
more  accurately  if  the  depth  were  decreased.   In  such  a  case 
the  problem  with  the  buoyancy  term  would  be  further  mini- 
mized. 

Further  studies  using  this  model  should  have  a  provi- 
sion for  including  vertical  compressibility.   Additionally, 
since  a  primary  source  of  the  energy  in  a  thunderstorm-type 
cloud  is  latent  heat,  application  to  cumulus  scale  convec- 
tive  activity  would  be  incomplete  without  the  inclusion  of 
moistureQ 
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APPENDIX  A 

EXTRAPOLATED    LIEBMANN   RELAXATION   SCHEME    FOR 
THREE   DIMENSIONAL   APPLICATION 

The   equation   to  be    solved  by  relaxation   is 

2 
V3     <t>      -   F   =    0  (A-l) 

where 

F  =    -  V  •      |(V,    •      V-)V,         +    fnf      + 


[^3    '      V3^3] 


3     |_v  3         3    3  I         °  eo    az 


Writing    (A-l)     in    finite   difference    form   and   expanding  yields 


<t>.  .  ~     .    ,     -20.      .    .     +  <f>.     0     .    , 
i+2,j,k  ri,  J,k  T  ^i-2,  j,k 

(2  Ax)2 

<*>i,j  +  2,k   -2<fti,j,k  +     fti,j-2,k 
(2  Ay)2 

<*>i,j,k+2   ~2  ^i,  j,k  +  ^i,  j,k-2 


-  F.     .    ,    =  0 


(2Az)2  i'j'k 


Rearranging   terms    leads    to 


*i,j,k=[       (2A^)2    (2    AZ)2    <*i+2,j,k  +    <*>i-2,j,k> 

+  (2Ax)2  (2Az)2  (0i/j  +  2/k  +  0i,j-2,k> 
+  (2Ax)2(2Ay)2  (<*>i/j/k+2  +  *i,JfM> 
-    (2  Ax)2(2Ay)2(2Az)2    Fi#JJ      A 
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where     A      is    the    fraction 

1 


2         (2Ay)2(2Az)2   +    (2Ax)2    (2Az)2   +    (2  Ax)  2  (2  Ay)  21 


old 

Adding  and  subtracting   0        produces 

i,  j,k 


new        old        r 
<*>.     .  .  =  <t>  .     .  _   +  A    (2Ay)2(2Az)2(  0.        +  0        ) 


+  (2Ax)2(2Az)2(«?».(.  +  2/k+0./._2/k) 


+  (2Ax)2(2Ay)2(«.ij/k+2+<?>i/jjk_2) 


i^-*  -  (2Ax)2(2Ay)2(2AZ)2Fi/j/J 


For  successive  over  relaxation  the  coefficient  (X      is 
included,  which  results  in  the  final  form 

new  old 

*i.j.k  =  (1-a>*i,j,k+  aA  [(2A^2(2Az>2(*i+2,j,k 

+  *i-2.J.k> 

+  (2Ax)2(2Az)2(0.ij  +  2/k+  0±/j_2/k) 

+  (2Ax)2(2Ay)2(^>iij/k+2  +0i/j/k_2) 

-  <2Ax)2(2Ay)2(2Az)2  Fi/JjkJ  (A-2) 
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Note    that   if    Ax  =    Ay  =    Az,    equation    (A-2)    reduces    to    the 
more    familiar 

new  old 

i-2, j,k 


new  o±a 


+     <*>i,j  +  2,k  +      M.,  j-2,k  +     C*)i/j,k+2   +     ^ijj,!^ 


-  <2*x>2  Fi.j,J  - 
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APPENDIX  B 
DEVELOPMENT  OF  THE  FORCING  FUNCTION 

In  order  to  maintain  consistent  finite  differencing 
between  the  prognostic  equations  and  the  forcing  function  of 
the  Liebmann  relaxation  equation,  a  procedure  which  was 
recommended  by  Dr.  Rc  T.  Williams  (personal  communication) 
is  used  to  develop  the  forcing  function.   The  process  is 
carried  out  in  6  steps  for  each  point  in  the  grid.   The 
6  steps  required  to  build   F  5  5  5,   the  forcing  function 
at  1=5,  j=5/  k=5,  are 

pi        ^(W)5.5.4 
5/5,5       2  Az 


2        -G(v)5/4/5       x 
5,5,5       2Ay         5,5,5 


-G(u)4  5  5       y 
5,5,5      2  Ax  5,5,5 


p4       __  G<">6,5,5      p3 

5,5,5      2  Ax  5,5,5 


p5         G(V)5,6,5      F4 

5,5,5      2  Ay         5,5,5 


F6       =  0(w)55,6   +   F5 

5,5,5       2Az  5,5,5 


where 
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G  (w)    =    -     £  (W)    + 


<30  2 


G(v)    =  -     £(v)    -   fu  +    y    v       v, 


and 


G(u)    =  -    £  (u)    +   fv  +   y    v       u. 
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APPENDIX  C 
BOUNDARY  CONDITIONS  AT  THE  WALLS 

For  the  plane  described  by  j=3,  the  equations  of  motion 
may  be  written 


1,  3,k  i, J/ k 


d  v  „  b<P 

--  £(v)i,3/k-r;        -  fuif3fk         (c-2) 


at  i,j,k      ay 

i/3,k  i,3,k 


^  -  -    tf(w)lf3f]c-22  +  g*i,3,k  .      Cc-3) 


dt  '    '         dz  6 

i,  3,k  i,3,k  o 


Forming  the  time  derivative  of  divergence  from  (C-l) ,  (C-2) , 
and  (C-3)  and  using  finite  differencing  yields 


du  du 


^  t  i+l,3,k               dt      i-l,3,k 
0   =  — 

2    Ax 


d  v  dv 


^t      i,4,k  ht      i,2,k 

2    Ay 

d  w  aw 


ht   i,j,k+l  ^t      j,j,k-l 

2    A  z 
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(C-4) 


m 


^educing  ^e 


condition 


Cc 


-5) 


vi,2.* 


-   -  vi,3,* 


^icn  i^lie! 


dt      if2,^ 


cc.4)    ^^ 


^3 


o  «= 


(C-6) 


of    CC-6)    V 


ields 


2*       ,   ^^  *      SCVU.3JJ. 

_.     £.^4-1,3.*      -     — -^       ~~2 

.  ^UiS^i  + --■-■- *** 

£  l*>  1,3.*!*. x)  ^ 


.    +   ui,3^     + 

2  &y 


2  A  z   a0 
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Applying    the    conditions 

ui,2,k  =   ui,3,k      '  (c"8> 

*i,3,k-*i,l.k=   "2Ay£ui,2,k      '  <C"9> 

and 

♦i.4.k-     *i,2,k=-2AyfUi,3,k      '  <C-10> 

produces 

^2*i,3,k=    Fi.3,k      • 

which  is  the  equation  to  be  solved  by  relaxation. 

Thus,  the  conditions  for  the  southern  boundary  must 
include  the  relationships  given  by  (C-5) ,  (C-8) ,  (C-9)  and 
(C-10) . 

A  similar  approach  produces  corresponding  conditions 
at  the  northern  boundary.   These  conditions  are 

vl,JM-l,k  =  "  vi,JM-2,k   *  (C-ll) 

ui,JM-l,k  =  ui,JM-2,k   '  (C-12) 

<t>.    _.  .  -  <f>.  n        =  -2Ayf  u-  _.  ,  .   ,     (C-13) 

l,JM,k    ^i,JM-2,k        x         lfJM-lfk   '     v     ' 


and 


*i.jM-i,k  -  *i,JM-3,k  ■  -2Ayf  ui,JM-2,k  •    <c-14) 
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For  the  plane  described  by  k=3,  the  same  process  can 
be  carried  out  to  develop  boundary  conditions  at  the  top  and 
bottom ,  which  are 


Wi,j#2  =  -Wi,j,3   <  (c"15> 


tu  ^  2Azg   . 


*i.j.4-  *i,J<2=i|f  •i.j.a  -      <c"17> 


and 


Wi,j,KM-l  ~  -wi,j,KM-2   '  (C-18> 


2  Azg 
^i,j,KM-"  *i,j,KM-2  "    0Q~  ^i#j,KM-l     (C-19) 


^4    I   ™   ,   -   ^. 


2  Azg 


i,j,KM-l    ^i,j,KM-3      0o    vi,j,KM-2   # 

(C-20) 
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